/* Copyright 2012,2013 Tobias Marschall
 * 
 * This file is part of CLEVER.
 * 
 * CLEVER is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * CLEVER is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with CLEVER.  If not, see <http://www.gnu.org/licenses/>.
 */

#include <stdio.h>
#include <execinfo.h>
#include <signal.h>
#include <stdlib.h>
#include <ctime>
#include <limits>

#include <iostream>

#include <boost/iostreams/device/file.hpp>
#include <boost/iostreams/filtering_stream.hpp>
#include <boost/iostreams/filter/gzip.hpp>
#include <boost/unordered_set.hpp>
#include <boost/unordered_map.hpp>
#include <boost/program_options.hpp>

#include <bamtools/api/BamWriter.h>

#include "VersionInfo.h"

#include "HistogramBasedDistribution.h"
#include "Histogram.h"
#include "NamedDnaSequence.h"
#include "FastaReader.h"
#include "FastqReader.h"
#include "SequenceUtils.h"
#include "GroupWiseBamReader.h"
#include "ShortDnaSequence.h"
#include "NfaMatcher.h"
#include "SplitAligner.h"
#include "BamHelper.h"
#include "split/Mapping.h"
#include "split/SplitMapping.h"
#include "split/WholeMapping.h"
#include "split/Read.h"
#include "VariationUtils.h"
#include "ThreadPool.h"
#include "MismatchWeightTrack.h"

using namespace std;
namespace po = boost::program_options;
namespace io = boost::iostreams;

typedef vector<NamedDnaSequence*> reference_list_t;

void usage(const char* name, const po::options_description& options_desc) {
	cerr << "Usage: " << name << " [options] <reference.fasta(.gz)> <unsplit.1.fastq(.gz)> <unsplit.2.fastq(.gz)>" << endl;
	cerr << "       " << name << " [options] --single-end <reference.fasta(.gz)> <unsplit.fastq(.gz)>" << endl;
	cerr << endl;
	cerr << "Reads BAM format of aligned split reads from stdin, joins them," << endl;
	cerr << "and outputs a alignments in BAM format to stdout." << endl;
	cerr << endl;
	cerr << "<insert-length-dist>: Insert length distribution as obtained by running" << endl;
	cerr << "                      insert-length-histogram on a BAM file." << endl;
	cerr << endl;
	cerr << options_desc << endl;
	exit(1);
}

/** Record on a complete read pair. */
typedef struct read_pair_record_t {
	string name;
	bool skip;
	Read read1;
	Read read2;
	read_pair_record_t(const string& name) : name(name), skip(false) {}
} read_pair_record_t;

typedef struct pairing_t {
	int aln1_idx;
	int aln2_idx;
	double p;
	int longest_indel;
	pairing_t(int aln1_idx, int aln2_idx, double p, int longest_indel) : aln1_idx(aln1_idx), aln2_idx(aln2_idx), p(p), longest_indel(longest_indel) {}
} pairing_t;

typedef struct {
	bool operator()(const pairing_t& p1, const pairing_t& p2) {
		if (p1.p != p2.p) return p1.p > p2.p;
		return p1.longest_indel < p2.longest_indel;
	}
} pairing_comparator_t;

/** Parameters and global data structures needed when processing a read. */
typedef struct parameters_t {
	const reference_list_t& references;
	const vector<string>& ref_names;
	const SplitAligner& aligner;
	HistogramBasedDistribution* insert_length_distribution;
	int max_mappings;
	int max_insert_size;
	int max_span;
	int anchor_length;
	int anchor_errors;
	int anchor_search_distance;
	int anchor_search_iterations;
	int max_anchors;
	int max_anchor_pairs;
	int anchor_length_increment;
	int max_anchor_length;
	bool verbose;
	bool simple_cigar;
	int phred_offset;
	bool store_mismatches;
	bool rightmost;
	bool interchromosomal;
	bool single_end;
	string read_group;
	parameters_t(const reference_list_t& references, const vector<string>& ref_names, const SplitAligner& aligner, HistogramBasedDistribution* insert_length_distribution, int max_mappings, int max_insert_size, int max_span, int anchor_length, int anchor_errors, int anchor_search_distance, int anchor_search_iterations, int max_anchors, int max_anchor_pairs, int anchor_length_increment, int max_anchor_length, bool verbose, bool simple_cigar, int phred_offset, bool store_mismatches, bool rightmost, bool interchromosomal, bool single_end, string read_group) : 
		references(references), ref_names(ref_names), aligner(aligner), insert_length_distribution(insert_length_distribution), max_mappings(max_mappings), max_insert_size(max_insert_size), max_span(max_span), anchor_length(anchor_length), anchor_errors(anchor_errors), anchor_search_distance(anchor_search_distance), anchor_search_iterations(anchor_search_iterations), max_anchors(max_anchors), max_anchor_pairs(max_anchor_pairs), anchor_length_increment(anchor_length_increment), max_anchor_length(max_anchor_length), verbose(verbose), simple_cigar(simple_cigar), phred_offset(phred_offset), store_mismatches(store_mismatches), rightmost(rightmost), interchromosomal(interchromosomal), single_end(single_end), read_group(read_group) {}
} parameters_t;

typedef struct weighted_variation_t {
	Variation variation;
	double weight;
	weighted_variation_t(Variation variation, double weight) : variation(variation), weight(weight) {}
} weighted_variation_t;

/** Values to be determined during processing of a read. */
typedef struct read_output_t {
	bool unique;
	int insert_size;
	vector<BamTools::BamAlignment> bam_aln1;
	vector<BamTools::BamAlignment> bam_aln2;
	ostringstream output;
	vector<weighted_variation_t> candidates;
	vector<MismatchWeightTrack::mismatch_weight_t> mismatches;
	bool valid_pair_found;
	bool skipped;
	read_output_t() : unique(false), insert_size(0), valid_pair_found(false), skipped(false) {}
} read_output_t;

/** Records the size of the largest indel in the given alignment in the given histogram. */
void record_indel_size(const BamTools::BamAlignment& alignment, Histogram* insertion_histogram, Histogram* deletion_histogram) {
	int largest_deletion = 0;
	int largest_insertion = 0;
	vector<BamTools::CigarOp>::const_iterator it = alignment.CigarData.begin();
	for (; it != alignment.CigarData.end(); ++it) {
		if (it->Type == 'D') {
			// don't record alignment if it contains multiple indels.
			if ((largest_deletion != 0) || (largest_insertion != 0)) return;
			largest_deletion = it->Length;
		}
		if (it->Type == 'I') {
			// don't record alignment if it contains multiple indels.
			if ((largest_deletion != 0) || (largest_insertion != 0)) return;
			largest_insertion = it->Length;
		}
	}
	if (insertion_histogram != 0) {
		insertion_histogram->add(largest_insertion);
	}
	if (deletion_histogram) {
		deletion_histogram->add(largest_deletion);
	}
}

void ensure_mapping_count(Read& read, read_output_t& output, const parameters_t& parameters, string& read_processing, int read_nr) {
	if ((int)read.getMappings().size() > parameters.max_mappings) {
		string filtering = "MR";
		int best_phred = read.getMappings()[0]->getPhredScore();
		read.filterByPhredScore(best_phred + parameters.phred_offset);
		if (parameters.verbose) output.output << "Read" << read_nr << ": too many mappings only retaining those with PHRED <= " << (best_phred + parameters.phred_offset) << endl;
		if ((int)read.getMappings().size() > parameters.max_mappings) {
			string filtering = "MB";
			read.filterByPhredScore(best_phred);
			if (parameters.verbose) output.output << "Read" << read_nr <<": still too many mappings only retaining optimal ones." << endl;
			if ((int)read.getMappings().size() > parameters.max_mappings) {
				string filtering = "MO";
				read.filterByPhredScore(best_phred - 1);
				if (parameters.verbose) output.output << "Read" << read_nr <<": still too many mappings, discarding all of them." << endl;
			}
		}
		if (read_processing.size() > 0) read_processing += ',';
		read_processing += filtering;
	}
}

void compute_read_posterior_distribution(Read& read, vector<double>& distribution) {
	vector<int> phreds(read.getMappings().size(), -1);
	int min_phred = numeric_limits<int>::max();
	for (size_t i=0; i<read.getMappings().size(); ++i) {
		phreds[i] = read.getMappings()[i]->getPhredScore();
		if (phreds[i] < min_phred) min_phred = phreds[i];
	}
	double prob_sum = 0.0;
	distribution.assign(read.getMappings().size(), 0.0);
	for (size_t i=0; i<read.getMappings().size(); ++i) {
		double p = pow(10.0, -((double)(phreds[i]-min_phred))/10.0);
		distribution[i] = p;
		prob_sum += p;
	}
	for (size_t i=0; i<read.getMappings().size(); ++i) {
		distribution[i] /= prob_sum;
	}
}

uint16_t probability_to_mapq(double probability) {
	if (probability <= 0.0) return (uint16_t)0;
	if (probability >= 1.0) return (uint16_t)100;
	return (uint16_t)round(min(-10.0*log1p(-probability)/log(10.0),100.0));
}

void finish_bam_record_paired_end(BamTools::BamAlignment& bam_aln, BamTools::BamAlignment* mate_aln, read_output_t& output, const parameters_t& parameters, bool is_primary, bool is_first_mate, bool mate_unmapped, double probability) {
	bam_aln.SetIsPrimaryAlignment(is_primary);
	bam_aln.SetIsPaired(true);
	bam_aln.MapQuality = probability_to_mapq(probability);
	bam_aln.SetIsFirstMate(is_first_mate);
	bam_aln.SetIsSecondMate(!is_first_mate);
	bam_aln.SetIsMapped(true);
	if (!bam_aln.AddTag("YP","f",(float)probability)) assert(false);
	if (mate_aln == 0) {
		bam_aln.SetIsMateMapped(!mate_unmapped);
	} else {
		assert(!mate_unmapped);
		bam_aln.SetIsProperPair(true);
		bam_aln.MateRefID = mate_aln->RefID;
		bam_aln.MatePosition = mate_aln->Position;
		bam_aln.SetIsMateMapped(true);
		bam_aln.SetIsMateReverseStrand(mate_aln->IsReverseStrand());
		if (bam_aln.Position < mate_aln->Position) {
			bam_aln.InsertSize = mate_aln->GetEndPosition(false, true) - bam_aln.Position + 1;
		} else {
			bam_aln.InsertSize = mate_aln->Position - bam_aln.GetEndPosition(false, true) - 1;
		}
	}
	// compute NM tag
	BamHelper::addNMTag(&bam_aln);
	// record mismatch weights, if requested
	if (parameters.store_mismatches) {
		MismatchWeightTrack::extractFromAlignment(bam_aln, probability, parameters.phred_offset, &output.mismatches);
	}
	// turn X and = CIGAR operations into Ms, if requested
	if (parameters.simple_cigar) {
		BamHelper::cigarReduceMismatches(&bam_aln.CigarData);
	}
	if (parameters.read_group.size() > 0) {
		if (!bam_aln.AddTag("RG", "Z", parameters.read_group)) assert(false);
	}
}

void finish_bam_record_single_end(BamTools::BamAlignment& bam_aln, read_output_t& output, const parameters_t& parameters, bool is_primary, double probability) {
	bam_aln.SetIsPrimaryAlignment(is_primary);
	bam_aln.SetIsPaired(false);
	bam_aln.MapQuality = probability_to_mapq(probability);
	bam_aln.SetIsFirstMate(false);
	bam_aln.SetIsSecondMate(false);
	bam_aln.SetIsMapped(true);
	// set is mate UNmapped to false (otherwise BAM does not validate with PICARD)
	bam_aln.SetIsMateMapped(true);
	// compute NM tag
	BamHelper::addNMTag(&bam_aln);
	if (!bam_aln.AddTag("YP","f",(float)probability)) assert(false);
	// record mismatch weights, if requested
	if (parameters.store_mismatches) {
		MismatchWeightTrack::extractFromAlignment(bam_aln, probability, parameters.phred_offset, &output.mismatches);
	}
	// turn X and = CIGAR operations into Ms, if requested
	if (parameters.simple_cigar) {
		BamHelper::cigarReduceMismatches(&bam_aln.CigarData);
	}
	if (parameters.read_group.size() > 0) {
		if (!bam_aln.AddTag("RG", "Z", parameters.read_group)) assert(false);
	}
}

void create_unmapped_bam_record(const Read& read, const string& name, BamTools::BamAlignment& bam_aln, BamTools::BamAlignment* mate_aln, bool is_first_mate, const parameters_t& parameters) {
	bam_aln.MapQuality = 0;
	bam_aln.SetIsPaired(true);
	bam_aln.SetIsFirstMate(is_first_mate);
	bam_aln.SetIsSecondMate(!is_first_mate);
	bam_aln.SetIsMapped(false);
	bam_aln.Name = name;
	bam_aln.QueryBases = read.getSequence().toString();
	bam_aln.Qualities = read.getSequence().qualityString();
	if (mate_aln != 0) {
		bam_aln.RefID = mate_aln->RefID;
		bam_aln.Position = mate_aln->Position;
		bam_aln.MateRefID = mate_aln->RefID;
		bam_aln.MatePosition = mate_aln->Position;
		bam_aln.SetIsMateMapped(true);
		bam_aln.SetIsMateReverseStrand(mate_aln->IsReverseStrand());
	}
	if (parameters.read_group.size() > 0) {
		if (!bam_aln.AddTag("RG", "Z", parameters.read_group)) assert(false);
	}
}

void create_unmapped_bam_record_single_end(const Read& read, const string& name, BamTools::BamAlignment& bam_aln, const parameters_t& parameters) {
	bam_aln.MapQuality = 0;
	bam_aln.SetIsPaired(false);
	bam_aln.SetIsFirstMate(false);
	bam_aln.SetIsSecondMate(false);
	bam_aln.SetIsMapped(false);
	// set is mate UNmapped to false (otherwise BAM does not validate with PICARD)
	bam_aln.SetIsMateMapped(true);
	bam_aln.Name = name;
	bam_aln.QueryBases = read.getSequence().toString();
	bam_aln.Qualities = read.getSequence().qualityString();
	if (parameters.read_group.size() > 0) {
		if (!bam_aln.AddTag("RG", "Z", parameters.read_group)) assert(false);
	}
}

void process_read_pair(read_pair_record_t& record, read_output_t& output, const parameters_t& parameters) {
	if (record.skip) {
		output.skipped = true;
		return;
	}
	if (parameters.verbose) {
		output.output << "============================================" << endl;
		output.output << "Processing read: " << record.name << endl;
		output.output << "--- INPUT" << endl;
		output.output << "Read1:" << endl;
		output.output << record.read1;
		output.output << "Read2:" << endl;
		output.output << record.read2;
	}

	// remarks about processing of reads to be put in YP tags
	// N: None / Normal processing
	// A<len>: increased anchor length to <len>
	// MR: only keeping alignments at most phred-offset worse than best
	// MB: only keeping best alignment(s)
	// MO: overflow in the number of best mappings
	string read1_processing = "";
	string read2_processing = "";
	
	int anchor_length = parameters.anchor_length;
	bool anchor_search_complete1 = false;
	bool anchor_search_complete2 = false;
	while (true) {
		if (parameters.verbose) output.output << "---- SEARCHING FOR LENGTH " << anchor_length << " ANCHORS"<< endl;
		Read::anchor_search_regions_t regions_already_searched1;
		Read::anchor_search_regions_t regions_already_searched2;
		for (int i=0; i<parameters.anchor_search_iterations; ++i) {
			if (parameters.verbose) output.output << "---- ANCHOR SEARCH ITERATION " << i << endl;
			auto_ptr<Read::anchor_search_regions_t> search_regions1(0);
			auto_ptr<Read::anchor_search_regions_t> search_regions2(0);
			if (!anchor_search_complete1) {
				search_regions1 = record.read1.determineAnchorSearchRegions(parameters.anchor_search_distance, anchor_length, &record.read2);
				search_regions1->subtract(regions_already_searched1);
				if (parameters.verbose) output.output << "Search regions for read1:" << endl;
				if (parameters.verbose) output.output << (*search_regions1);
			}
			if (!anchor_search_complete2) {
				search_regions2 = record.read2.determineAnchorSearchRegions(parameters.anchor_search_distance, anchor_length, &record.read1);
				search_regions2->subtract(regions_already_searched2);
				if (parameters.verbose) output.output << "Search regions for read2:" << endl;
				if (parameters.verbose) output.output << (*search_regions2);
			}
			if (!anchor_search_complete1) {
				record.read1.findNewAnchors(*search_regions1, parameters.references, anchor_length, parameters.anchor_errors);
				regions_already_searched1.add(*search_regions1);
			}
			if (!anchor_search_complete2) {
				record.read2.findNewAnchors(*search_regions2, parameters.references, anchor_length, parameters.anchor_errors);
				regions_already_searched2.add(*search_regions2);
			}
		}
		if (!anchor_search_complete1) {
			if (parameters.verbose) output.output << "Anchors for read1: left: " << record.read1.leftAnchors().size() << ", right: " << record.read1.rightAnchors().size() << ", pairs: " << (record.read1.leftAnchors().size() * record.read1.rightAnchors().size()) << endl;
			if ((record.read1.leftAnchors().size() <= parameters.max_anchors) && (record.read1.rightAnchors().size() <= parameters.max_anchors) &&
				(record.read1.leftAnchors().size() * record.read1.rightAnchors().size() <= parameters.max_anchor_pairs)
			) {
				if (parameters.verbose) output.output << "Anchors for read1: search complete" << endl;
				anchor_search_complete1 = true;
			} else {
				record.read1.resetAnchorList();
			}
		}
		if (!anchor_search_complete2) {
			if (parameters.verbose) output.output << "Anchors for read2: left: " << record.read2.leftAnchors().size() << ", right: " << record.read2.rightAnchors().size() << ", pairs: " << (record.read2.leftAnchors().size() * record.read2.rightAnchors().size()) << endl;
			if ((record.read2.leftAnchors().size() <= parameters.max_anchors) && (record.read2.rightAnchors().size() <= parameters.max_anchors) &&
				(record.read2.leftAnchors().size() * record.read2.rightAnchors().size() <= parameters.max_anchor_pairs)
			) {
				if (parameters.verbose) output.output << "Anchors for read2: search complete" << endl;
				anchor_search_complete2 = true;
			} else {
				record.read2.resetAnchorList();
			}
		}
		if (anchor_search_complete1 && anchor_search_complete2) {
			break;
		} else {
			assert(parameters.anchor_length_increment > 0);
			anchor_length += parameters.anchor_length_increment;
			if (anchor_length > parameters.max_anchor_length) {
				// Read ends still having too many anchors will be ignored: we treat them as if no anchor was found
				if (!anchor_search_complete1) record.read1.clearAnchorList();
				if (!anchor_search_complete2) record.read2.clearAnchorList();
				break;
			} else {
				if (parameters.verbose) output.output << "Increasing anchor length to " << anchor_length << endl;
			}
		}
	}
	
	if (parameters.verbose) output.output << "--- COMPUTE ALIGNMENTS FOR READ 1" << endl;
	record.read1.computeAlignments(parameters.references, parameters.aligner, parameters.max_span, parameters.interchromosomal, parameters.max_anchor_pairs);
	if (parameters.verbose) output.output << "--- COMPUTE ALIGNMENTS FOR READ 2" << endl;
	record.read2.computeAlignments(parameters.references, parameters.aligner, parameters.max_span, parameters.interchromosomal, parameters.max_anchor_pairs);

	// sort mappings by phred score
	record.read1.sortMappings();
	record.read2.sortMappings();
	bool unique = (record.read1.getMappings().size() == 1) && (record.read2.getMappings().size() == 1);
	if (parameters.verbose) {
		output.output << "--- RESULTING READ RECORD" << endl;
		output.output << "Read1:" << endl;
		output.output << record.read1;
		output.output << "Read2:" << endl;
		output.output << record.read2;
	}

	ensure_mapping_count(record.read1, output, parameters, read1_processing, 1);
	ensure_mapping_count(record.read2, output, parameters, read2_processing, 2);

	// compute pairings for alignments of both read ends
	vector<pairing_t> pairings;
	double prob_sum = 0.0;
	for (size_t i=0; i<record.read1.getMappings().size(); ++i) {
		for (size_t j=0; j<record.read2.getMappings().size(); ++j) {
			const Mapping& mapping1 = *record.read1.getMappings()[i];
			const Mapping& mapping2 = *record.read2.getMappings()[j];
			if (mapping1.getRefId() != mapping2.getRefId()) continue;
			if (mapping1.isReverse() == mapping2.isReverse()) continue;
			if (mapping1.getStartPosition() == mapping2.getStartPosition()) continue;
			const Mapping& left = (mapping1.getStartPosition()<mapping2.getStartPosition())?mapping1:mapping2;
			const Mapping& right = (mapping1.getStartPosition()<mapping2.getStartPosition())?mapping2:mapping1;
			if (!right.isReverse()) continue;
			int insert_length = right.getStartPosition() - left.getEndPosition() - 1;
			if (insert_length > parameters.max_insert_size) continue;
			double p = left.getProbability() * right.getProbability();
			if (parameters.insert_length_distribution != 0) {
				p *= parameters.insert_length_distribution->probability(insert_length);
			}
			if (p == 0.0) continue;
			prob_sum += p;
			pairings.push_back(pairing_t(i,j,p, max(left.longestIndel(),right.longestIndel())));
		}
	}
	
	// Probability of each (end) alignment being correct,
	vector<double> probs1(record.read1.getMappings().size(), 0.0);
	vector<double> probs2(record.read2.getMappings().size(), 0.0);
	// Index of mate of each alignment (-1 if unpaired)
	vector<int> mates1(record.read1.getMappings().size(), -1);
	vector<int> mates2(record.read2.getMappings().size(), -1);

	if (pairings.size() == 0) {
		if (parameters.verbose) output.output << "No valid pair found!" << endl;
		// if no pairing is available, compute probabilites for each read individually
		compute_read_posterior_distribution(record.read1, probs1);
		compute_read_posterior_distribution(record.read2, probs2);
	} else {
		output.valid_pair_found = true;
		// sort pairings by probability
		sort(pairings.begin(), pairings.end(), pairing_comparator_t());
		// normalize probabilites, i.e. compute posteriors
		for (size_t i=0; i<pairings.size(); ++i) {
			pairings[i].p /= prob_sum;
			if (parameters.verbose) output.output << " Valid pairing: " << pairings[i].aln1_idx << ", " << pairings[i].aln2_idx << ", " << pairings[i].p << endl;
		}
		// Greedily pair all reads with the most probable unpaired mate...
		for (size_t i=0; i<pairings.size(); ++i) {
			probs1[pairings[i].aln1_idx] += pairings[i].p;
			probs2[pairings[i].aln2_idx] += pairings[i].p;
			// Skip pairing if one of the reads is already paired
			if ((mates1[pairings[i].aln1_idx] != -1) || (mates2[pairings[i].aln2_idx] != -1)) continue;
			// store pair
			mates1[pairings[i].aln1_idx] = pairings[i].aln2_idx;
			mates2[pairings[i].aln2_idx] = pairings[i].aln1_idx;
		}
	}
	// if pair is uniquely mapped, record some statistics
	if (unique) {
		// if insert size distribution is to be estimated, record this read
		const Mapping& mapping1 = *record.read1.getMappings()[0];
		const Mapping& mapping2 = *record.read2.getMappings()[0];
		const Mapping& left = (mapping1.getStartPosition()<mapping2.getStartPosition())?mapping1:mapping2;
		const Mapping& right = (mapping1.getStartPosition()<mapping2.getStartPosition())?mapping2:mapping1;
		output.unique = true;
		output.insert_size = right.getStartPosition() - left.getEndPosition() - 1;
		if (parameters.verbose) output.output << "Read pair maps uniquely --> recording insert size " << output.insert_size << endl;
	}

	// Generate all BAM records (to be modified later)
	for (size_t i=0; i<record.read1.getMappings().size(); ++i) {
		output.bam_aln1.push_back(record.read1.getMappings()[i]->getBamAlignment(record.read1, record.name));
		// Store all putative variations
		auto_ptr<vector<Variation> > variations = record.read1.getMappings()[i]->getPutativeVariations(parameters.references, parameters.rightmost);
		if (variations.get() != 0) {
			for (vector<Variation>::const_iterator it = variations->begin(); it !=variations->end(); ++it) {
				output.candidates.push_back(weighted_variation_t(*it, probs1[i]));
			}
		}
	}
	for (size_t i=0; i<record.read2.getMappings().size(); ++i) {
		output.bam_aln2.push_back(record.read2.getMappings()[i]->getBamAlignment(record.read2, record.name));
		// Store all putative variations
		auto_ptr<vector<Variation> > variations = record.read2.getMappings()[i]->getPutativeVariations(parameters.references, parameters.rightmost);
		if (variations.get() != 0) {
			for (vector<Variation>::const_iterator it = variations->begin(); it !=variations->end(); ++it) {
				output.candidates.push_back(weighted_variation_t(*it, probs2[i]));
			}
		}
	}
	// Finish alignment(s) for first read
	if (record.read1.getMappings().size() == 0) {
		// Read is unmapped, create "unmapped" record
		output.bam_aln1.push_back(BamTools::BamAlignment());
		create_unmapped_bam_record(record.read1, record.name, output.bam_aln1[0], (output.bam_aln2.size()==0)?0:(&output.bam_aln2[0]), true, parameters);
	} else {
		bool mate_unmapped = record.read2.getMappings().size() == 0;
		for (size_t i=0; i<record.read1.getMappings().size(); ++i) {
			BamTools::BamAlignment& bam_aln = output.bam_aln1[i];
			bool is_primary = false;
			if (pairings.size() > 0) {
				is_primary = (int)i == pairings[0].aln1_idx;
			} else {
				is_primary = i == 0;
			}
			finish_bam_record_paired_end(bam_aln, (mates1[i]==-1)?0:&output.bam_aln2[mates1[i]], output, parameters, is_primary, true, mate_unmapped, probs1[i]);
		}
	}
	// Finish alignment(s) for second read
	if (record.read2.getMappings().size() == 0) {
		// Read is unmapped, create "unmapped" record
		output.bam_aln2.push_back(BamTools::BamAlignment());
		create_unmapped_bam_record(record.read2, record.name, output.bam_aln2[0], (output.bam_aln1.size()==0)?0:(&output.bam_aln1[0]), false, parameters);
	} else {
		bool mate_unmapped = record.read1.getMappings().size() == 0;
		for (size_t i=0; i<record.read2.getMappings().size(); ++i) {
			BamTools::BamAlignment& bam_aln = output.bam_aln2[i];
			bool is_primary = false;
			if (pairings.size() > 0) {
				is_primary = (int)i == pairings[0].aln2_idx;
			} else {
				is_primary = i == 0;
			}
			finish_bam_record_paired_end(bam_aln, (mates2[i]==-1)?0:&output.bam_aln1[mates2[i]], output, parameters, is_primary, false, mate_unmapped, probs2[i]);
		}
	}
}

void process_single_read(read_pair_record_t& record, read_output_t& output, const parameters_t& parameters) {
	if (record.skip) {
		output.skipped = true;
		return;
	}
	if (parameters.verbose) {
		output.output << "============================================" << endl;
		output.output << "Processing read: " << record.name << endl;
		output.output << "--- INPUT" << endl;
		output.output << "Read1:" << endl;
		output.output << record.read1;
		output.output << "Read2:" << endl;
		output.output << record.read2;
	}

	// remarks about processing of reads to be put in YP tags
	// N: None / Normal processing
	// A<len>: increased anchor length to <len>
	// MR: only keeping alignments at most phred-offset worse than best
	// MB: only keeping best alignment(s)
	// MO: overflow in the number of best mappings
	string read1_processing = "";
	
	int anchor_length = parameters.anchor_length;
	bool anchor_search_complete1 = false;
	while (true) {
		if (parameters.verbose) output.output << "---- SEARCHING FOR LENGTH " << anchor_length << " ANCHORS"<< endl;
		Read::anchor_search_regions_t regions_already_searched1;
		for (int i=0; i<parameters.anchor_search_iterations; ++i) {
			if (parameters.verbose) output.output << "---- ANCHOR SEARCH ITERATION " << i << endl;
			auto_ptr<Read::anchor_search_regions_t> search_regions1(0);
			if (!anchor_search_complete1) {
				search_regions1 = record.read1.determineAnchorSearchRegions(parameters.anchor_search_distance, anchor_length, 0);
				search_regions1->subtract(regions_already_searched1);
				if (parameters.verbose) output.output << "Search regions for read1:" << endl;
				if (parameters.verbose) output.output << (*search_regions1);
				record.read1.findNewAnchors(*search_regions1, parameters.references, anchor_length, parameters.anchor_errors);
				regions_already_searched1.add(*search_regions1);
			}
		}
		if (!anchor_search_complete1) {
			if (parameters.verbose) output.output << "Anchors for read1: left: " << record.read1.leftAnchors().size() << ", right: " << record.read1.rightAnchors().size() << ", pairs: " << (record.read1.leftAnchors().size() * record.read1.rightAnchors().size()) << endl;
			if ((record.read1.leftAnchors().size() <= parameters.max_anchors) && (record.read1.rightAnchors().size() <= parameters.max_anchors) &&
				(record.read1.leftAnchors().size() * record.read1.rightAnchors().size() <= parameters.max_anchor_pairs)
			) {
				if (parameters.verbose) output.output << "Anchors for read1: search complete" << endl;
				anchor_search_complete1 = true;
			} else {
				record.read1.resetAnchorList();
			}
		}
		if (anchor_search_complete1) {
			break;
		} else {
			assert(parameters.anchor_length_increment > 0);
			anchor_length += parameters.anchor_length_increment;
			if (anchor_length > parameters.max_anchor_length) {
				// Read ends still having too many anchors will be ignored: we treat them as if no anchor was found
				if (!anchor_search_complete1) record.read1.clearAnchorList();
				break;
			} else {
				if (parameters.verbose) output.output << "Increasing anchor length to " << anchor_length << endl;
			}
		}
	}
	
	if (parameters.verbose) output.output << "--- COMPUTE ALIGNMENTS FOR READ 1" << endl;
	record.read1.computeAlignments(parameters.references, parameters.aligner, parameters.max_span, parameters.interchromosomal, parameters.max_anchor_pairs);

	// sort mappings by phred score
	record.read1.sortMappings();
	bool unique = (record.read1.getMappings().size() == 1);
	if (parameters.verbose) {
		output.output << "--- RESULTING READ RECORD" << endl;
		output.output << "Read1:" << endl;
		output.output << record.read1;
	}

	ensure_mapping_count(record.read1, output, parameters, read1_processing, 1);
	
	// Probability of each (end) alignment being correct,
	vector<double> probs1(record.read1.getMappings().size(), 0.0);
	compute_read_posterior_distribution(record.read1, probs1);
	// if readr is uniquely mapped, record some statistics
	if (unique) {
		output.unique = true;
	}

	// Generate all BAM records
	if (record.read1.getMappings().size() == 0) {
		// Read is unmapped, create "unmapped" record
		output.bam_aln1.push_back(BamTools::BamAlignment());
		create_unmapped_bam_record_single_end(record.read1, record.name, output.bam_aln1[0], parameters);
	} else {
		output.valid_pair_found = true;
		for (size_t i=0; i<record.read1.getMappings().size(); ++i) {
			output.bam_aln1.push_back(record.read1.getMappings()[i]->getBamAlignment(record.read1, record.name));
			finish_bam_record_single_end(output.bam_aln1[output.bam_aln1.size()-1], output, parameters, i==0, probs1[i]);
			// Store all putative variations
			auto_ptr<vector<Variation> > variations = record.read1.getMappings()[i]->getPutativeVariations(parameters.references, parameters.rightmost);
			if (variations.get() != 0) {
				for (vector<Variation>::const_iterator it = variations->begin(); it !=variations->end(); ++it) {
					output.candidates.push_back(weighted_variation_t(*it, probs1[i]));
				}
			}
		}
	}
}

void printBadXTagsMessage(GroupWiseBamReader* bam_reader, const string& filename, bool skip_non_xa) {
	cerr << "Reading from \"" << filename << "\":" << endl;
	cerr << bam_reader->getInconsistentXTagsMessage() << endl;
	cerr << "  This can have several reasons, for example:" << endl;
	cerr << "  1) You use BWA and have forgotten to run xa2multi.pl (distributed with BWA)." << endl;
	cerr << "  2) The used read aligner does not write X0 and X1 tags correctly." << endl;
	if (!skip_non_xa) {
		cerr << "  3) If you use BWA alignments, remember to use option -X." << endl;
	}
	cerr << endl;
	cerr << "To ignore this inconsistency, use option -I." << endl;
}

bool isSoftClipped(const BamTools::BamAlignment& a) {
	if (!a.IsMapped()) return false;
	for (size_t i=0; i<a.CigarData.size(); ++i) {
		if (a.CigarData[i].Type == 'S') return true;
	}
	return false;
}

bool checkReferences(const GroupWiseBamReader& bam_reader, const FastaReader::reference_map_t& ref_map, reference_list_t* ref_list) {
	const BamTools::RefVector& ref_data = bam_reader.getReferenceData();
	assert(ref_list->size() == 0);
	for (size_t i=0; i<ref_data.size(); ++i) {
		FastaReader::reference_map_t::const_iterator ref_it = ref_map.find(ref_data[i].RefName);
		if (ref_it == ref_map.end()) {
			cerr << "Error: Reference sequence for \"" << ref_data[i].RefName << "\" not found." << endl;
			return false;
		}
		ref_list->push_back(ref_it->second);
	}
	return true;
}

void read_unsplit_paired_end(FastqReader& reader1, FastqReader& reader2, auto_ptr<FastqReader::fastq_record_t>* target1, auto_ptr<FastqReader::fastq_record_t>* target2) {
	if (!reader1.hasNext()) throw runtime_error("Unexpected end of FASTQ input 1.");
	if (!reader2.hasNext()) throw runtime_error("Unexpected end of FASTQ input 2.");
	*target1 = reader1.getNext();
	*target2 = reader2.getNext();
	if ((*target1)->name.compare((*target2)->name) != 0) {
		ostringstream oss;
		oss << "FASTQ input files out of sync. Read name " << (*target1)->name << " != " << (*target2)->name;
		throw runtime_error(oss.str());
	}
}

void read_unsplit_single_end(FastqReader& reader, auto_ptr<FastqReader::fastq_record_t>* target) {
	if (!reader.hasNext()) throw runtime_error("Unexpected end of FASTQ input.");
	*target = reader.getNext();
}

typedef struct work_package_t {
	read_pair_record_t* record;
	const parameters_t& parameters;
	read_output_t output;

	work_package_t(read_pair_record_t* record, const parameters_t& parameters) : record(record), parameters(parameters) {}
	void run() {
		if (parameters.single_end) {
			process_single_read(*record, output, parameters);
		} else {
			process_read_pair(*record, output, parameters);
		}
	}
} work_package_t;

typedef boost::unordered_map<Variation,double> variation_weight_map_t;

typedef struct output_writer_t {
	BamTools::BamWriter& bam_writer;
	variation_weight_map_t* variation_candidates;
	Histogram* insert_size_estimator;
	Histogram* insertion_histogram;
	Histogram* deletion_histogram;
	MismatchWeightTrack* mismatches_tracks;
	bool output_secondary;
	long long totally_processed;
	long long valid_pair_found_count;
	long long skipped;

	output_writer_t(BamTools::BamWriter& bam_writer, variation_weight_map_t* variation_candidates, Histogram* insert_size_estimator, Histogram* insertion_histogram, Histogram* deletion_histogram, MismatchWeightTrack* mismatches_tracks, bool output_secondary) : bam_writer(bam_writer), variation_candidates(variation_candidates), insert_size_estimator(insert_size_estimator), insertion_histogram(insertion_histogram), deletion_histogram(deletion_histogram), mismatches_tracks(mismatches_tracks), output_secondary(output_secondary), totally_processed(0), valid_pair_found_count(0), skipped(0) {}
	
	void write(auto_ptr<work_package_t> work) {
		if (variation_candidates != 0) {
			vector<weighted_variation_t>::const_iterator it = work->output.candidates.begin();
			for (; it != work->output.candidates.end(); ++it) {
				variation_weight_map_t::iterator map_it = variation_candidates->find(it->variation);
				if (map_it == variation_candidates->end()) {
					(*variation_candidates)[it->variation] = it->weight;
				} else {
					map_it->second += it->weight;
				}
			}
		}
		if (mismatches_tracks != 0) {
			mismatches_tracks->addAll(work->output.mismatches);
		}
		totally_processed += 1;
		if (work->output.skipped) skipped += 1;
		if (work->output.valid_pair_found) valid_pair_found_count += 1;

		if (work->output.unique) {
			if (insert_size_estimator != 0) insert_size_estimator->add(work->output.insert_size);
			// if indel size distribution is to be estimated, record this read
			if ((deletion_histogram != 0) || (insertion_histogram != 0)) {
				assert(work->output.bam_aln1.size() == 1);
				assert(work->output.bam_aln2.size() == 1);
				record_indel_size(work->output.bam_aln1[0], insertion_histogram, deletion_histogram);
				record_indel_size(work->output.bam_aln2[0], insertion_histogram, deletion_histogram);
			}
		}
		for (size_t i=0; i<work->output.bam_aln1.size(); ++i) {
			const BamTools::BamAlignment& a = work->output.bam_aln1[i];
			if (output_secondary || a.IsPrimaryAlignment()) {
				bam_writer.SaveAlignment(a);
			}
		}
		if (!work->parameters.single_end) {
			for (size_t i=0; i<work->output.bam_aln2.size(); ++i) {
				const BamTools::BamAlignment& a = work->output.bam_aln2[i];
				if (output_secondary || a.IsPrimaryAlignment()) {
					bam_writer.SaveAlignment(a);
				}
			}
		}
		cerr << work->output.output.str();
		delete work->record;
	}
} output_writer_t;

void signal_handler(int sig) {
  void *array[20];
  size_t size;
  size = backtrace(array, 20);
  fprintf(stderr, "Error: signal %d:\n", sig);
  backtrace_symbols_fd(array, size, 2);
  exit(1);
}

int main(int argc, char* argv[]) {
	signal(SIGSEGV, signal_handler);
	signal(SIGABRT, signal_handler);
	VersionInfo::checkAndPrintVersion("laser-core", cerr);
	string commandline = VersionInfo::commandline(argc, argv);

	// PARAMETERS
	bool skip_non_xa = false;
	bool ignore_wrong_x_tags = false;
	int max_mappings;
	int cost_threshold;
	int indel_cost;
	bool affine_gapcosts = false;
	int gap_extend_cost;
	int split_cost;
	int interchrom_split_cost;
	int secondary_aln_phred_diff;
	int phred_offset;
	int anchor_length;
	int anchor_errors;
	int anchor_search_distance;
	int anchor_search_iterations;
	int max_span;
	int max_anchors;
	int max_anchor_pairs;
	int max_insert_size;
	int anchor_length_increment;
	int max_anchor_length;
	int min_anchor_length;
	int threads;
	bool output_secondary = false;
	string putative_variations_filename = "";
	string insert_length_histogram_input_filename = "";
	string insert_length_histogram_output_filename = "";
	string insertion_length_histogram_filename = "";
	string deletion_length_histogram_filename = "";
	bool verbose = false;
	bool simple_cigar = false;
	string snp_filename = "";
	double snp_weight_cutoff;
	double indel_weight_cutoff;
	int max_input_aln = -1;
	bool rightmost = false;
	bool interchromosomal = false;
	bool single_end = false;
	bool allow_soft_clipping = false;
	int soft_clip_open_cost;
	int soft_clip_extend_cost;
	string read_group = "";
	string read_group_library = "";
	string read_group_sample = "";

	po::options_description options_desc("Allowed options");
	options_desc.add_options()
		("verbose,v", po::value<bool>(&verbose)->zero_tokens(), "Be (very) verbose. CAUTION: produces a lot of output to stderr.")
		("single-end", po::value<bool>(&single_end)->zero_tokens(), "Process single end reads (instead of paired end).")
		("skip_non_xa,X", po::value<bool>(&skip_non_xa)->zero_tokens(), "Skip reads for which other alignments exist (i.e. X0+X1>1), but no XA tag is present. Turn on when using BWA.")
		("max_input_aln", po::value<int>(&max_input_aln)->default_value(-1), "Maximum number of input anchor alignments. If more are given, skip this anchor (but not the whole read).")
		("ignore_wrong_X_tags,I", po::value<bool>(&ignore_wrong_x_tags)->zero_tokens(), "Do not abort when wrongly set X0/X1 tags are encountered.")
		("phred_offset,p", po::value<int>(&phred_offset)->default_value(33), "Value to subtract from ASCII code to get the PHRED quality.")
		("max_mapping,m", po::value<int>(&max_mappings)->default_value(10), "Maximum number of mappings for each read in a pair. If more mappings are found, the read is skipped.")
		("cost_threshold,c", po::value<int>(&cost_threshold)->default_value(115), "Maximum PHRED sum allowed for (split) read alignments.")
		("indel_cost,i", po::value<int>(&indel_cost)->default_value(30), "PHRED like cost of an indel.")
		("affine_gapcosts", po::value<bool>(&affine_gapcosts)->zero_tokens(), "Use affine gapcosts (instead of linear gapcosts). If enabled, gap open costs are given by --indel_cost and gap extend costs are given by --gap_extend_cost, i.e. a gap of length k costs indel_cost+k*gap_extend_cost.")
		("gap_extend_cost", po::value<int>(&gap_extend_cost)->default_value(10), "Gap extend cost. Only effective in connection with --affine_gapcosts.")
		("split_cost,s", po::value<int>(&split_cost)->default_value(25), "PHRED like cost of splitting an alignment.")
		("interchromosomal", po::value<bool>(&interchromosomal)->zero_tokens(), "Search for interchromosomal split reads.")
		("interchrom_split_cost", po::value<int>(&interchrom_split_cost)->default_value(60), "PHRED like cost of an inter-chromosomal split.")
		("secondary_aln_phred_diff,d", po::value<int>(&secondary_aln_phred_diff)->default_value(29), "Report/evaluate split alignments with a difference of at most this value to the best split alignment.")
		("min_anchor_length,a", po::value<int>(&min_anchor_length)->default_value(8), "Minimal number of nucleotides on each sides of a split.")
		("max_span,n", po::value<int>(&max_span)->default_value(1000), "Maximal allowed span (i.e. length on reference) of a split read alignment.")
		("max_anchors", po::value<int>(&max_anchors)->default_value(100), "Maximal allowed anchors per read suffix/prefix.")
		("max_anchor_pairs", po::value<int>(&max_anchor_pairs)->default_value(500), "Maximal allowed number of anchor pairings.")
		("anchor_length_increment", po::value<int>(&anchor_length_increment)->default_value(4), "Number of characters to increase anchor length by if too many anchors are found.")
		("max_anchor_length", po::value<int>(&max_anchor_length)->default_value(30), "Maximal anchor length (if reached it will not be increased any further.")
		("max_insert,N", po::value<int>(&max_insert_size)->default_value(50000), "Maximal allowed internal segment size when pairing alignments for a read pair.")
		("output_secondary,S", po::value<bool>(&output_secondary)->zero_tokens(), "Also output secondary alignments.")
		("anchor_search_length,A", po::value<int>(&anchor_length)->default_value(20), "Length of anchor used for searching.")
		("anchor_errors,e", po::value<int>(&anchor_errors)->default_value(2), "Allowed errors when searching for anchors.")
		("anchor_distance,o", po::value<int>(&anchor_search_distance)->default_value(600), "Regions size for searching for anchors.")
		("anchor_search_iter,t", po::value<int>(&anchor_search_iterations)->default_value(3), "Anchor search iterations.")
		("putative_variations,P", po::value<string>(&putative_variations_filename)->default_value(""), "Filename to write a list of all variation candidates to. All candidates are written, no matter how weak the evidence for them to be true.")
		("input_length_dist_in,l", po::value<string>(&insert_length_histogram_input_filename)->default_value(""), "File with known internal segment length histogram to be used to score alignments.")
		("input_length_dist_out,L", po::value<string>(&insert_length_histogram_output_filename)->default_value(""), "File to write internal segment length histogram for uniquely mappable reads to.")
		("insertion_length_dist,R", po::value<string>(&insertion_length_histogram_filename)->default_value(""), "File to write insertion length histogram for uniquely mappable reads to.")
		("deletion_length_dist,D", po::value<string>(&deletion_length_histogram_filename)->default_value(""), "File to write deletion length histogram for uniquely mappable reads to.")
		("threads,T", po::value<int>(&threads)->default_value(0), "Number of threads (default: 0 = strictly single-threaded).")
		("m_in_cigar,M", po::value<bool>(&simple_cigar)->zero_tokens(), "Use M for matches and mismatches in CIGAR strings (instead of '=' and 'X').")
		("snp", po::value<string>(&snp_filename)->default_value(""), "Filename to store putative SNP positions to.")
		("snp_weight_cutoff,w", po::value<double>(&snp_weight_cutoff)->default_value(3.0), "Weight cutoff for SNPs to be written to filename given as parameter --snp.")
		("indel_weight_cutoff", po::value<double>(&indel_weight_cutoff)->default_value(3.0), "Weight cutoff for indels to be written to filename given as parameter -P.")
		("rightmost", po::value<bool>(&rightmost)->zero_tokens(), "Report rightmost deletions/insertions (default: leftmost).")
		("soft_clip", po::value<bool>(&allow_soft_clipping)->zero_tokens(), "Allow reads to be soft-clipped.")
		("soft_clip_open_cost", po::value<int>(&soft_clip_open_cost)->default_value(35), "Cost for soft clipping a read.")
		("soft_clip_extend_cost", po::value<int>(&soft_clip_extend_cost)->default_value(3), "Cost for \"extending\" a soft clip; i.e., softclipping k characters from a read will cost soft_clip_open_cost+k*soft_clip_extend_cost.")
		("read_group", po::value<string>(&read_group)->default_value(""), "Read group for all reads.")
		("read_group_library", po::value<string>(&read_group_library)->default_value(""), "Library information for read group header.")
		("read_group_sample", po::value<string>(&read_group_sample)->default_value(""), "Sample information for read group header.")
	;

	for (int i=1; i<argc; ++i) {
		string arg(argv[i]);
		if (arg.compare("--single-end") == 0) single_end = true;
	}
	if (isatty(fileno(stdin)) || (single_end && (argc<3)) || (!single_end && (argc<4))) {
		usage(argv[0], options_desc);
	}

	string reference_filename = "";
	string unsplit1_filename = "";
	string unsplit2_filename = "";

	if (single_end) {
		reference_filename = argv[argc-2];
		unsplit1_filename = argv[argc-1];
		argc -= 2;
	} else {
		reference_filename = argv[argc-3];
		unsplit1_filename = argv[argc-2];
		unsplit2_filename = argv[argc-1];
		argc -= 3;
	}

	po::variables_map options;
	try {
		po::store(po::parse_command_line(argc, argv, options_desc), options);
		po::notify(options);
	} catch(exception& e) {
		cerr << "error: " << e.what() << "\n";
		return 1;
	}

	if (single_end) {
		if (insert_length_histogram_input_filename.size() > 0) {
			cerr << "Error: options --single_end and -l cannot be combined." << endl;
			return 1;
		}
		if (insert_length_histogram_output_filename.size() > 0) {
			cerr << "Error: options --single_end and -L cannot be combined." << endl;
			return 1;
		}
	}

	cerr << "Commandline: " << commandline << endl;

	SplitAligner aligner(cost_threshold, (affine_gapcosts?SplitAligner::AFFINE:SplitAligner::LINEAR), indel_cost, gap_extend_cost, split_cost, interchrom_split_cost, min_anchor_length, secondary_aln_phred_diff, phred_offset, false, rightmost, allow_soft_clipping, soft_clip_open_cost, soft_clip_extend_cost);
	
	auto_ptr<FastaReader::reference_map_t> reference_map(0);
	try {
		reference_map = FastaReader::parseFromFile(reference_filename);
		cerr << "Read " << reference_map->size() << " reference sequences." << endl;
	}catch(exception& e) {
		cerr << "Error: " << e.what() << "\n";
		return 1;
	}

	ifstream* unsplit1_istream = 0;
	io::filtering_istream* in1 = 0;
	FastqReader* fastq_reader1 = 0;
	ifstream* unsplit2_istream = 0;
	io::filtering_istream* in2 = 0;
	FastqReader* fastq_reader2 = 0;

	unsplit1_istream = new ifstream(unsplit1_filename.c_str());
	if (unsplit1_istream->fail()) {
		cerr << "Error opening file \"" + unsplit1_filename + "\"." << endl;
		return 1;
	}
	in1 = new io::filtering_istream();
	if (unsplit1_filename.substr(unsplit1_filename.size()-3,3).compare(".gz") == 0) {
		in1->push(io::gzip_decompressor());
	}
	in1->push(*unsplit1_istream);
	fastq_reader1 = new FastqReader(*in1, true);

	if (!single_end) {
		unsplit2_istream = new ifstream(unsplit2_filename.c_str());
		if (unsplit2_istream->fail()) {
			cerr << "Error opening file \"" + unsplit2_filename + "\"." << endl;
			return 1;
		}
		in2 = new io::filtering_istream();
		if (unsplit2_filename.substr(unsplit2_filename.size()-3,3).compare(".gz") == 0) {
			in2->push(io::gzip_decompressor());
		}
		in2->push(*unsplit2_istream);
		fastq_reader2 = new FastqReader(*in2, true);
		
	}

	MismatchWeightTrack* mismatches_tracks = 0;
	ofstream* snps_outfile = 0;
	bool store_mismatches = snp_filename.size() > 0;
	if (store_mismatches) {
		mismatches_tracks = new MismatchWeightTrack();
		snps_outfile = new ofstream(snp_filename.c_str());
		if (snps_outfile->fail()) {
			cerr << "Error: could not open " << snp_filename << " for writing." << endl;
			return 1;
		}
	}

	// If output of all putative variations is requested (option -P), then
	// initialize set and respective output stream
	variation_weight_map_t* variation_candidates = 0;
	ofstream* variation_candidates_file = 0;
	if (putative_variations_filename.size() > 0) {
		variation_candidates = new variation_weight_map_t();
		variation_candidates_file = new ofstream(putative_variations_filename.c_str());
		if (variation_candidates_file->fail()) {
			cerr << "Error: could not open " << variation_candidates_file << " for writing." << endl;
			return 1;
		}
	}

	// prepare estimation and output of insert size distribution if requested
	Histogram* insert_length_histogram = 0;
	ofstream* insert_length_distribution_outfile = 0;
	if (insert_length_histogram_output_filename.size() > 0) {
		insert_length_distribution_outfile = new ofstream(insert_length_histogram_output_filename.c_str());
		if (insert_length_distribution_outfile->fail()) {
			cerr << "Error: could not open " << insert_length_distribution_outfile << " for writing." << endl;
			return 1;
		}
		insert_length_histogram = new Histogram();
	}

	// prepare estimation and output of insertion lengths if requested
	Histogram* insertion_length_histogram = 0;
	ofstream* insertion_length_distribution_outfile = 0;
	if (insertion_length_histogram_filename.size() > 0) {
		insertion_length_distribution_outfile = new ofstream(insertion_length_histogram_filename.c_str());
		if (insertion_length_distribution_outfile->fail()) {
			cerr << "Error: could not open " << insertion_length_histogram_filename << " for writing." << endl;
			return 1;
		}
		insertion_length_histogram = new Histogram();
	}

	// prepare estimation and output of deletion lengths if requested
	Histogram* deletion_length_histogram = 0;
	ofstream* deletion_length_distribution_outfile = 0;
	if (deletion_length_histogram_filename.size() > 0) {
		deletion_length_distribution_outfile = new ofstream(deletion_length_histogram_filename.c_str());
		if (deletion_length_distribution_outfile->fail()) {
			cerr << "Error: could not open " << deletion_length_histogram_filename << " for writing." << endl;
			return 1;
		}
		deletion_length_histogram = new Histogram();
		// put one pseudocount on the longest possible deletion
		deletion_length_histogram->add(max_span);
	}
	clock_t clock_start = clock();
	GroupWiseBamReader* splitread_bam_reader = 0;
	HistogramBasedDistribution* insert_length_distribution = 0;
	auto_ptr<FastqReader::fastq_record_t> whole_read1(0);
	auto_ptr<FastqReader::fastq_record_t> whole_read2(0);
	output_writer_t* output_writer = 0;
	try {
		splitread_bam_reader = new GroupWiseBamReader("/dev/stdin", false);
//		splitread_bam_reader = new GroupWiseBamReader("input.bam", false);
		splitread_bam_reader->enableProgressMessages(cerr, 200000);
		reference_list_t reference_list;
		if (!checkReferences(*splitread_bam_reader,*reference_map,&reference_list)) return 1;
		// prepare vector with names of reference sequences
		vector<string> ref_names;
		for (size_t i = 0; i<reference_list.size(); ++i) {
			ref_names.push_back(reference_list[i]->getName());
		}
		// initialize BamWriter to write alignments to stdout
		BamTools::BamWriter* bam_writer = 0;
		bam_writer = new BamTools::BamWriter();
		BamTools::SamHeader bam_output_header(splitread_bam_reader->getHeader());
		BamTools::SamProgram laser_core_program("laser-core");
		laser_core_program.CommandLine = commandline;
		laser_core_program.Version = VersionInfo::version();
		laser_core_program.Name = "laser-core";
		bam_output_header.Programs.Add(laser_core_program);
		if (read_group.size() > 0) {
			BamTools::SamReadGroup rg(read_group);
			if (read_group_library.size() > 0) rg.Library = read_group_library;
			if (read_group_sample.size() > 0) rg.Sample = read_group_sample;
			bam_output_header.ReadGroups.Add(rg);
		}
		bam_writer->Open("/dev/stdout", bam_output_header, splitread_bam_reader->getReferenceData());
		//bam_writer->Open("/dev/null", splitread_bam_reader->getHeader(), splitread_bam_reader->getReferenceData());
		// load insert size distribution from file if requested
		if (insert_length_histogram_input_filename.size() > 0) {
			insert_length_distribution = new HistogramBasedDistribution(insert_length_histogram_input_filename);
		}
		// Initialize thread pool
		typedef ThreadPool<work_package_t,output_writer_t> thread_pool_t;
		output_writer = new output_writer_t(*bam_writer, variation_candidates, insert_length_histogram, insertion_length_histogram, deletion_length_histogram, mismatches_tracks, output_secondary);
		thread_pool_t* thread_pool = new thread_pool_t(threads, 100, threads, *output_writer);
		// Assemble set of parameters to be used to process reads
		parameters_t parameters(reference_list, ref_names, aligner, insert_length_distribution, max_mappings, max_insert_size, max_span, anchor_length, anchor_errors, anchor_search_distance, anchor_search_iterations, max_anchors, max_anchor_pairs, anchor_length_increment, max_anchor_length, verbose, simple_cigar, phred_offset, store_mismatches, rightmost, interchromosomal, single_end, read_group);
		auto_ptr<read_pair_record_t> current_pair(new read_pair_record_t(""));
		while ( splitread_bam_reader->hasNext() ) {
			splitread_bam_reader->advance();
			if (splitread_bam_reader->getReadName().size() < 4) {
				cerr << "Illegal read name: \"" << splitread_bam_reader->getReadName() << "\"" << endl;
				return 1;
			}
			if (!ignore_wrong_x_tags && splitread_bam_reader->hasInconsistentXTags()) {
				printBadXTagsMessage(splitread_bam_reader, "/dev/stdin", skip_non_xa);
				return 1;
			}
			const string& n = splitread_bam_reader->getReadName();
			string original_name = n.substr(0, n.size()-3);
			if (original_name.compare(current_pair->name) != 0) {
				if (current_pair->name.size()>0) {
					auto_ptr<work_package_t> wp(new work_package_t(current_pair.release(), parameters));
					thread_pool->addTask(wp);
// 					process_read_pair(*current_pair, reference_list, ref_names, aligner, insert_length_distribution, max_span, max_mappings, max_insert_size, bam_writer, output_secondary, variation_candidates, insert_length_histogram, insertion_length_histogram, deletion_length_histogram, verbose);
				}
				current_pair = auto_ptr<read_pair_record_t>(new read_pair_record_t(original_name));
				if (single_end) {
					read_unsplit_single_end(*fastq_reader1, &whole_read1);
				} else {
					read_unsplit_paired_end(*fastq_reader1, *fastq_reader2, &whole_read1, &whole_read2);
				}
				if (whole_read1->name.compare(original_name) != 0) {
					cerr << "Error: input out of sync. Read name \"" << original_name << "\" from stdin, but \"" << whole_read1->name << "\" from FASTQ input." << endl;
					return 1;
				}
				if (!current_pair->skip) {
					current_pair->read1.addSequenceAndQualities(whole_read1->sequence.toString(),whole_read1->sequence.qualityString());
					if (!single_end) {
						current_pair->read2.addSequenceAndQualities(whole_read2->sequence.toString(),whole_read2->sequence.qualityString());
					}
				}
			}
			string name_suffix = n.substr(n.size()-3, 3);
			// cerr << "  ... read name suffix: " << name_suffix << endl;
			if ((name_suffix[0]!='_') || ((name_suffix[1]!='1') && (name_suffix[1]!='2')) || ((name_suffix[2]!='R') && (name_suffix[2]!='L'))){
				cerr << "Illegal read name: \"" << splitread_bam_reader->getReadName() << "\"" << endl;
				return 1;
			}
			//const vector<ReadRecord::anchor_t>& anchors = (name_suffix[1] == '1')?current_pair->read1.anchors:current_pair->read2.anchors;
			bool first_read = name_suffix[1] == '1';
			bool left = name_suffix[2] == 'L';
			if (single_end && !first_read) {
				cerr << "Error: encountered read with suffix " << name_suffix << ", which is illegal for single end reads" << endl;
				return 1;
			}
			Read& read = first_read?current_pair->read1:current_pair->read2;
			typedef vector<BamTools::BamAlignment*>::const_iterator bam_aln_it_t;
			bool skip = false;
			if (skip_non_xa && splitread_bam_reader->hasMultipleMappingsButNoXA()) {
				if (verbose) cerr << "Warning: ignoring alignments for " << splitread_bam_reader->getReadName() << ": no XA tag" << endl;
				skip = true;
			}
			if ((max_input_aln != -1) && (splitread_bam_reader->getAlignments().size() > max_input_aln)) {
				if (verbose) cerr << "Warning: ignoring alignments for " << splitread_bam_reader->getReadName() << ": found " << splitread_bam_reader->getAlignments().size() << " alignments, but " << max_input_aln << " is maximum (option --max_input_aln)" << endl;
				skip = true;
			}
			if (!skip) {
				for (bam_aln_it_t it = splitread_bam_reader->getAlignments().begin(); it != splitread_bam_reader->getAlignments().end(); ++it) {
					if (!(*it)->IsMapped()) continue;
					if ((*it)->RefID < 0) continue;
					if ((*it)->GetEndPosition() - 1 >= (int)reference_list[(*it)->RefID]->size()) {
						if (verbose) cerr << "Warning: Ignoring alignment that streches over reference end for " << (*it)->Name << endl;
						continue;
					}
					if (isSoftClipped(**it)) {
						if (verbose) cerr << "Warning: Ignoring soft-clipped split read alignment for \"" << n << "\"." << endl;
						continue;
					}
					if (left) {
						if ((*it)->IsReverseStrand()) {
							read.addRightAnchor((*it)->RefID, (*it)->GetEndPosition()-1, true, (*it)->Length, false);
						} else {
							read.addLeftAnchor((*it)->RefID, (*it)->Position, false, (*it)->Length, false);
						}
					} else {
						if ((*it)->IsReverseStrand()) {
							read.addLeftAnchor((*it)->RefID, (*it)->Position, true, (*it)->Length, false);
						} else {
							read.addRightAnchor((*it)->RefID, (*it)->GetEndPosition()-1, false, (*it)->Length, false);
						}
					}
				}
			}
		}
		if (current_pair->name.size()>0) {
			auto_ptr<work_package_t> wp(new work_package_t(current_pair.release(), parameters));
			thread_pool->addTask(wp);
// 			process_read_pair(*current_pair, reference_list, ref_names, aligner, insert_length_distribution, max_span, max_mappings, max_insert_size, bam_writer, output_secondary, variation_candidates, insert_length_histogram, insertion_length_histogram, deletion_length_histogram, verbose);
		}
		delete thread_pool;
		if (bam_writer != 0) {
			bam_writer->Close();
			delete bam_writer;
		}
		if (fastq_reader1 != 0) delete fastq_reader1;
		if (in1 != 0) delete in1;
		if (unsplit1_istream != 0) delete unsplit1_istream;
		if (fastq_reader2 != 0) delete fastq_reader2;
		if (in2 != 0) delete in2;
		if (unsplit2_istream != 0) delete unsplit2_istream;
	} catch(exception& e) {
		cerr << "Error: " << e.what() << "\n";
		return 1;
	}
	// Write putative variations to file if requested
	if (variation_candidates != 0) {
		vector<Variation> variation_list;
		for (variation_weight_map_t::const_iterator it=variation_candidates->begin(); it!=variation_candidates->end(); ++it) {
			variation_list.push_back(it->first);
		}
		sort(variation_list.begin(), variation_list.end(), VariationUtils::variation_position_sort_t());
		vector<Variation>::const_iterator it = variation_list.begin();
		long long indels_above_threshold = 0;
		for (; it != variation_list.end(); ++it) {
			variation_weight_map_t::const_iterator map_it = variation_candidates->find(*it);
			assert(map_it != variation_candidates->end());
			if (map_it->second < indel_weight_cutoff) continue;
			indels_above_threshold += 1;
			(*variation_candidates_file) << (*it) << " " << map_it->second << endl;
		}
		cerr << indels_above_threshold << " of " << variation_list.size() << " putative indels have an expected support >= " << indel_weight_cutoff << " and are written to \"" << putative_variations_filename << "\"." << endl;
		variation_candidates_file->close();
	}
	double cpu_time = (double)(clock() - clock_start) / CLOCKS_PER_SEC;
	cerr << "CPU time used in alignment phase: " << cpu_time << endl;
	// Write insert size distribution to file
	if (insert_length_histogram != 0) {
		if (insert_length_histogram->size() < 1000000) {
			cerr << "Warning: less than 1000000 uniquely mappable reads found (" << insert_length_histogram->size() << "). Estimate of " << endl;
			cerr << "         insert size distribution written to \"" << insert_length_histogram_output_filename <<  "\" might be inaccurate." << endl;
		}
		(*insert_length_distribution_outfile) << (*insert_length_histogram->toDistribution(20));
		insert_length_distribution_outfile->close();
	}
	// Write insertion size distribution to file
	if (insertion_length_histogram != 0) {
		(*insertion_length_distribution_outfile) << (*insertion_length_histogram->toDistribution(100));
		insertion_length_distribution_outfile->close();
	}
	// Write deletion size distribution to file
	if (deletion_length_histogram != 0) {
		(*deletion_length_distribution_outfile) << (*deletion_length_histogram->toDistribution(100));
		deletion_length_distribution_outfile->close();
	}
	// Write SNPs to filename
	if (snps_outfile != 0) {
		auto_ptr<vector<MismatchWeightTrack::snp_t> > snps = mismatches_tracks->getSnpCandidates(snp_weight_cutoff);
		const BamTools::RefVector& ref_data = splitread_bam_reader->getReferenceData();
		for (size_t i=0; i<snps->size(); ++i) {
			const MismatchWeightTrack::snp_t& s = snps->at(i);
			(*snps_outfile) << ref_data[s.chromosome_id].RefName << ' ' << (s.position+1) << ' ' << s.weight << endl;
		}
		cerr << "Found " << snps->size() << " putative SNPs with an expected support >= " << snp_weight_cutoff << " and wrote them to \"" << snp_filename << "\"." << endl;
	}
	if (single_end) {
		cerr << "Processed a total of " << output_writer->totally_processed << " single end reads." << endl;
	} else {
		cerr << "Processed a total of " << output_writer->totally_processed << " read pairs." << endl;
	}
	if (output_writer->skipped > 0) {
		cerr << "Skipped " << output_writer->skipped << " reads." << endl;
	}
	if (single_end) {
		cerr << "Found valid alignments for " << output_writer->valid_pair_found_count << " single end reads." << endl;
	} else {
		cerr << "Found valid alignment pairs for " << output_writer->valid_pair_found_count << " read pairs." << endl;
	}
	delete output_writer;
	FastaReader::reference_map_t::const_iterator ref_it = reference_map->begin();
	for (; ref_it != reference_map->end(); ++ref_it) {
		delete ref_it->second;
	}
	if (insert_length_histogram != 0) delete insert_length_histogram;
	if (insert_length_distribution_outfile != 0) delete insert_length_distribution_outfile;
	if (splitread_bam_reader != 0) delete splitread_bam_reader;
	if (insert_length_distribution != 0) delete insert_length_distribution;
	if (variation_candidates != 0) delete variation_candidates;
	if (variation_candidates_file != 0) delete variation_candidates_file;
	return 0;
}
